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SUMMARY 


An analytical technique is presented for approximating unsteady aerody- 
namic forces in the time domain. The order of elements of a matrix Pade approx- 
imation is postulated, and the resulting polynomial coefficients are determined 
through a combination of least-squares estimates for the numerator coefficients 
and a constrained gradient search for denominator coefficients which insures 
stable approximating functions. The number of differential equations required 
to represent the aerodynamic forces to a given accuracy tends to be smaller than 
that employed in certain existing techniques where the denominator coefficients 
are chosen a priori. The resulting Pade approximation allows the aircraft equa- 
tions of motion to be formulated in terms of linear ordinary differential equa- 
tions, a form amenable to the application of optimal control theory methodology. 
The technique is applied to an aeroelastic, cantilevered, semispan wing whose 
motion is expressed in terms of five elastic modes. A good estimate is obtained 
to the generalized unsteady aerodynamic forces on the imaginary axis for ele- 
ments of the matrix Pade approximation having fourth-order numerator and second- 
order denominator polynomials. 


INTRODUCTION 

In the past two decades, control-system designers have seen a shift in the 
mathematical techniques employed for the analysis and synthesis of dynamic sys- 
tems. The emphasis has shifted from the so-called "classical” or frequency 
domain methods developed prior to 1960 by investigators such as Nyquist and 
Bode , which have played an important role with respect to single input/single 
output systems, to the more "modern" state-space approach. This shift in empha- 
sis can be attributed to the control designer not only requiring control but 
also needing or desiring to "optimize" the performance of the control system. 

A difficulty in using modern control theory for the design of systems to control 
aeroelastic behavior is the requirement of transforming the unsteady aerodynamic 
forces, normally provided in the frequency domain, into the time domain. Fur- 
thermore, it is highly desirable to use a transformation that results in a 
dynamic system that is of "reduced order," since the cost of the design 
increases rapidly with increases in the number of differential equations used 
to model the system. Several methods that have been presented in the literature 
to perform this transformation are outlined. 

A "Power Series Expansion Method" is presented by Weiss, Tseng, and Morino 
in reference 1 . This method assumes that the unsteady aerodynamic forces can 
be represented by a power- series expansion in frequency. An advantage of this 
method is that if the unsteady aerodynamic forces can be represented with two 
or less terms, then no additional state equations are added to the set of equa- 
tions used to represent the motion of the aircraft. However, as more terms are 
added in the expansion, higher order derivatives of the aircraft state variables 
are introduced. Since unsteady aerodynamic forces are characteristically pro- 
portional to the delayed aircraft state, the use of these higher order terms to 



approximate the unsteady aerodynamic forces generally will not be as satisfac- 
tory as a more appropriate approximation function for the same number of added 
state equations. 

Edwards has proposed, in reference 2, the use of a "Rational Model" to 
approximate the unsteady aerodynamic forces which add no additional state equa- 
tions to the mathematical model. This method has been applied only to simple 
examples because of the lack of a production computer code for the generation 
of the aerodynamic forces in the Laplace domain. 

The method most commonly employed is referred to as the "Least-Squares 
Method." This method has been used for a number of different aerodynamic con- 
figurations in references 3, 4, and 5. The method derives its name from the 
method used to solve a set of simultaneous equations for the coefficients of 
an assumed aerodynamic model. This assumed aerodynamic model consists of a 
rational polynomial in the frequency domain. In order to make the problem 
linear, the coefficients of the denominator polynomial are assumed known. A 
least-squares estimator is then employed to solve for the coefficients of the 
numerator polynomial. These coefficients minimize the mean-square error 
between the predicted aerodynamic forces and the aerodynamic force data being 
fit at a fixed set of frequencies. A disadvantage of this method, as it is 
applied, is the requirement of the user to have a priori knowledge of the denom- 
inator polynomial in the assumed aerodynamic model. These parameters, which 
are usually arbitrarily chosen to be within the dynamic range of the natural 
frequencies of the aircraft, determine the best accuracy that can be achieved 
with a given model order. The order of the aerodynamic model is directly pro- 
portional to the number of first-order constant coefficient differential equa- 
tions generated by the algorithm. For the control-systems designer to increase 
the accuracy of the approximation, he must either spend a considerable amount 
of time and effort adjusting the parameters of the denominator polynomial or 
increase the system order and pay the cost of additional state equations in the 
design cycle. 

Another method that has received attention in the literature is the "Matrix 
Pade Method" as described in references 2, 6, and 7. This method has been 
applied to a number of simple problems in references 2 and 7. However, when 
this method is applied to problems of increasing difficulty where there are more 
than two modes and/or when the predicted aerodynamics are not precise, the 
resulting approximation function can be unstable. This characteristic of the 
"Matrix Pade Method" was observed in reference 6. In the present paper, a func- 
tion is considered stable if a bounded input (i.e., wing motion) always produces 
a bounded output (i.e., distributed load on the airfoil). This must not be con- 
fused with a stable aerodynamic approximation function which, when coupled with 
the dynamics of the structure, results in a closed- loop system that has a closed- 
loop instability (i.e., flutter). 

The method described in the present paper is an extension of the "Matrix 
Pade Method" and the "Least-Squares Method." A constrained gradient method is 
used to find the denominator coefficients of the approximation while constrain- 
ing the calculated approximation to be stable, and a linear least-squares method 
is used to solve for the numerator coefficients. By using the method described 
herein, the control designer need only adjust the order of the approximation 
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until a satisfactory compromise is reached between the good accuracy obtained 
when a large number of states are used and the high cost of performing modern 
control- theory design on large systems, 

A brief development of the equations of motion for an elastic airplane are 
presented. The technique of approximating the unsteady aerodynamic forces is 
then described. This technique is applied to an aeroelastic, cantilevered, 
semispan-wing flutter model. Results obtained are compared with those found 
using the "Least-Squares Method" in terms of both the fit of the oscillatory 
aerodynamic forces and the predicted behavior of critical system characteristic 
roots in the vicinity of flutter. The method developed in the present paper 
utilizes a lower order approximation for the aerodynamic forces than that 
employed in the "Least-Squares Method" of reference 3. The present method tends 
in general to require a lower order of approximation than the "Least-Squares 
Method" because of the inclusion of the denominator polynomials as additional 
variables in the fit of the oscillatory aerdynamic forces. This results in a 
smaller number of differential equations for the design model. 


A 

[ Ai ] 



a i 

[B] 


tc] 

c r 

c/2 

[d] 

[D m ] 


SYMBOLS 

matrix of coefficients used in solution of aerodynamic approximation 
(see appendix) 

coefficient matrices of aerodynamic approximation of equation (21 ) 
amplitude of generalized aerodynamic force Qij (t) 

matrix of coefficients of first-order differential equations of 
complete system 

coefficients of numerator polynomial of a Pade approximation 
matrix of coefficients (see matrix [c] ) 

vector of constant coefficients used in solution of aerodynamic 
approximation (see appendix) 

coefficients of denominator polynomial of a Pade approximation 

matrix of coefficients for system equations such that 
[c] X(t) = [b] X(t) 

wing root chord (see fig. 2) 

reference length 

generalized damping matrix 

coefficient matrices of aerodynamic approximation of equation (21) 
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e natural base 

P[x(t) ] ,F -1 [x( ja>) ] Fourier transform pair defined by equation (9) 

[i] identity matrix 

J measure of the error between data and approximation 

Jj function minimized by equation (20) 

Jj measure of the error between data and approximation of jth vibration 

mode 

j 

K order of numerator polynomial of a Pade approximation 

[k] generalized stiffness matrix 

k reduced frequency, wc/2V 

specific reduced frequency 
M free-stream Mach number 

[m] generalized mass matrix 

Mj generalized mass for ith vibration mode 

m(x,y) mass per unit area at point x,y 

N order of unsteady aerodynamic Pade approximation and order of denomi- 

nator polynomial in Pade approximation 

Ni order of unsteady aerodynamic Pade approximation for ith vibration 

node 

Nfc number of reduced frequencies used in cost function 

n number of vibration modes 

[p m ] matrix of coefficients of [P(j(D) ] (see eq. (17)) 

[P m# j] jth column vector of [p m ] 

p ij (j k &) ith, jth element of [p] evaluated at frequency k& 

[p ( jto) ] , [p ] matrix of numerator polynomials of unsteady aerodynamic Pade 
approximation 

AP(x,y,t) time varying total pressure distribution 
4 
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APj (x,y, t) time varying pressure distribution due to impulse of jth mode 

[Q (t) ] matrix of functions such that ith, jth element is unit impulse 

response of ith generalized force due to an impulse of jth mode 

Qi(t) generalized aerodynamic force in time domain for ith mode 

Qij (t) ith, jth element of [Q(t) ] 

[Q(jw)] Fourier transform of [Q(t)] 

CQ ( joj) ] unsteady aerodynamic Pade approximation 

Qij ( jkj,) ith, jth element of [Q(jo))] evaluated at frequency k£ 

[Qij(jfy,)] ith, jth element of [Q( jw) ] evaluated at frequency k£ 
q(t) vector of generalized coordinate in time domain 

qj (t) ith generalized coordinate in time domain 

q(jw) Fourier transform of q(t) 

[RgJ matrix of coefficients of R(jui) (see eq. (16)) 

[r ( joi) ] matrix [R(jw)] with leading coefficient equal to [i] 

Rj(jty,) element of [R(ja))] evaluated at frequency k& 

[r( jco) ] , [ r] matrix of denominator polynomials of unsteady aerodynamic forces 
Pade approximation 

r& f j jth diagonal element of [r^] 

s semispan length of wing in example 

T(jk) general Pade approximation 

t time 

V free-stream velocity 

Xj solution vector used in calculation of aerodynamic approximation (see 

appendix) 

X (t) system state vector 

x(t)i state variable vector 

(x,y) scalar product 

z(x,y,t) vertical displacement of body at point x,y 
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P 

4>i 

<t>ij 

0) 

II II 

Dots over 


viscous damping coefficient for ith mode 
reference density of fluid 

nondimens ional mode shape used to define generalized coordinate 
qi(t) 

phase angle of Qij (t) 
frequency of oscillation 

natural frequency of vibration of ith mode 
norm, \j( , ) 

symbols denote derivatives with respect to time. 


EQUATIONS OF MOTION FOR AEROELASTIC MODEL 

The equations of motion for an elastic airplane in straight and level 
flight can be formulated by using Lagrange's equations of motion with orthogonal 
modes. By the method of separating variables, the motion on the wing is assumed 
to be the product of a position function and a time function. It is also 
assumed that the product of these functions can be represented with sufficient 
accuracy by a finite series so that the vertical displacement at the point x,y 
becomes 


z(x,y,t) 


l 

1=1 


qi(t) <j>i (x,y) 


(1) 


where <j>i(x,y) are the nondimensional mode shapes used to represent the system, 
n is the number of modes included, and qi(t) represents the generalized 
coordinate with dimension of length, the particular unit depending upon the 
system of units being used. If all of the structural damping present in the 
aircraft is assumed to be viscous in nature, then the equations of motion can 
be written as 


Mi qi<t) + 2C i M i io n . (t) + u>n.Mi q ± (t) = 


- - PV 2 Qi(t) 


( 2 ) 


where co n ^ is the natural frequency of vibration of the ith mode, is the 

viscous damping coefficient for the ith mode, and Mi is the generalized mass 
of the ith mode defined by 
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(3) 




m(x,y) <j>i(x,y) dx dy 


,y>] 


in which 

51 n ^ 

S 


s 

m(x,y) is the mass density of the body at the point x,y and 
dy is the integral of the function over the area of the body denoted 


by S. The term on the right-hand side of equation (2) represents the general- 
ized aerodynamic force applied to the ith mode and is defined by 


- pv 2 Qj_ (t) = ^ pv 2 jy 

s 


AP(x,y,t) 

1 2 
- pv 2 
2 


4>i(x,y) 


dx dy 


(4) 


where the term PV^/2 is the dynamic pressure of the free stream, Ap(x,y,t) 
is the time varying perturbation in the lifting pressure distribution at the 
point x,y, and ( t ) is the generalized aerodynamic force normalized by the 

dynamic pressure. The lifting pressure distribution will be assumed to be 
representable by the following equation: 


AP (x, y, t) 


Ap(x,y,0) 


t 

* 1 4 " j (x f y,t-T) 


qj(T) 


dx 


(5) 


where Apj(x,y,t) is the time history of the pressure at the point x,y due 
to the jth mode being displaced by the unit impulse function at t = 0. Substi- 
tuting this expression for Ap(x,y,t) into equation (4) with the initial pertur- 
bation in the pressure distribution equal to zero and factoring out the dynamic 
pressure, Qi(t) can be expressed as 


Qi(t) 


11 < 

-H 

j=i 


Qij(t-T) qj(T) dT 


(6) 


where Qij (t) is defined by 
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mm 


Qij(t) 



APj (x f y,t) 



4>i(x,y) 


dx dy 


(7) 


Equation (2) can be rewritten in matrix form as 


1 pt 

[M] q(t) + [D] 5(t) + [K] q(t) + - pV 2 \ [ Q ( t-T ) ] q(T) dT = 0 (8) 

2 Jq 


where [M] is a diagonal matrix with the ith diagonal element being defined by 
equation (3) , [D] is a diagonal matrix with the ith diagonal element being 2 
2 Ci M iWnj/ [K] is a diagonal matrix with the ith diagonal element being 

[ Q (t) ] is a matrix with the ith,jth element being defined by equation (7), q(t) 
is a vector with the ith element being the variable q^(t), and £ is a null 
matrix. 

Introducing the Fourier transform pair as 


_ 00 

x(jw) = F[x (t) ] = r x(t)e“3 wt dt 

•-A-oo 


i r°° 

:(t) = F -1 [x(jw)] = — \ xjjcoje^t dw 

2lT J_oo J 


(9) 


and applying the transform to equation (8) yields 


( jo>) 2 [M] + (jw) [D] + [K] + - pV 2 [Q ( jw) ] > q(jco) = 0 

2 


(10) 


where the following Fourier transform properties have been used: 
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(jw) 2 q(joj) = F[q(t) ] 
(jco) q( j(D) = F[q(t) ] 
q(jw) = F[q(t) ] 


> v 


( 11 ) 


fQ ( jto) ] q ( jco) 



fQ(t-T)] q(T) dT 


-i J 


Equation (10) represents a variation of the classical flutter equation, and is 
valid only for stable aircraft motion. 

The ith,jth element of the matrix [Q (jco) ] represents the Fourier trans- 
form of the ith generalized unsteady aerodynamic force, normalized by the 
dynamic pressure, due to a unit impulse of the jth mode. Alternatively, it can 
be interpreted as the contribution to the ith generalized force due to a steady 
oscillation in the jth mode at a frequency of w, which can be expressed in the 
time domain as 


Qij (t) — Aij sin (cot + 4*ij) 


( 12 ) 


where A^j 
of [Q(t)j 


and <J>ij are the amplitude and phase angle of the ith, jth element 
or [Q(]0))]. 


DESCRIPTION OF TECHNIQUE 

In this section an approximation is developed for the unsteady aerodynamic 
forces in the frequency domain and, by inversion, in the time domain. The 
approximation gives an accurate representation of the aerodynamic forces for 
oscillatory motion provided the assumptions of small perturbation, inviscid 
flow are valid. The primary sources of errors in the approximation are the 
availability of only a limited frequency band of aerodynamic force data due to 
oscillatory motion and the limited order of the matrix Pade approximation 
employed to fit these data. Quantification of the accuracy of the aerodynamic 
forces for arbitrary motion requires additional analytical and experimental 
research and is beyond the scope of this paper. 

Theoretical three-dimensional oscillatory aerodynamic forces are normally 
calculated at a specified Mach number with the body oscillating at a number of 
reduced frequencies. The reduced frequency k is a nondimens ional number that 
represents the number of radians through which the body oscillates per reference 
length c/2 the body travels through the fluid, or 


1 U)C 
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where a) is the dimensional circular frequency of oscillation and V is the 
velocity of the body through the fluid. Vepa, in reference 7, and Edwards, in 
reference 2, suggest that the unsteady aerodynamic forces may be approximated 
by a rational function of polynomials in the complex frequency domain (i.e., 
a Pade approximation in the reduced frequency domain) . Pade approximations are 
classified by the order of the polynomials of the numerator and denominator. 

In the present paper, the notation [K,N] is used to represent a Pade approxi- 
mation with a numerator polynomial of order K over a denominator polynomial 
of order N, or 


T(jk) 


I 

-i =n 


(jk) 1 ai 


I 


(jk)^bg, 


0 4) 


where T(jk) is the approximation of the function, a£ are the coefficients of 
the numerator polynomial, b& are the coefficients of the denominator (the coef- 
ficient bfj is set to unity) , and k is the independent variable of the approx- 
imation. Vepa, in reference 7, and Edwards, in reference 2, have suggested the 
use of an [n+1 ,n] Pade approximation based on the high-frequency asymptotic 
behavior of unsteady aerodynamics as predicted by piston theory. Since the 
[n+1 ,n] approximation is a special case of the [n+2,n] Pade approximation and 
since it is desired to approximate the unsteady aerodynamic forces only over a 
specified frequency range (thereby ignoring the high-frequency behavior) , an 
[n+2,n] Pade approximation was employed for the development of the technique 
described herein. This is also consistent with the order of the approximation 
used in references 3 to 5 for the "Least-Squares Method." 

The objective is to find a set of stable [n+2,n] Pade approximations 
which best fit the matrix [Q(jw)] of equation (10) at a discrete set of fre- 
quencies. The parameter N represents the number of terms used to represent 
the lag in the development of the circulation (for the two-dimensional case they 
would approximate the Theodorsen circulation function) ; henceforth, N is 
referred to as the order of the unsteady aerodynamic Pade approximation. Also, 
the term "Pade approximation" refers to the unsteady aerodynamic Pade approxima- 
tion described in the present paper. A necessary and sufficient condition for 
stability of the Pade approximation is that the roots of the denominator poly- 
nomial of the Pade approximation have negative real parts. Furthermore, it is 
appropriate, from the number of differential equations generated by the algo- 
rithm and the amount of computational resources required by the algorithm, to 
require that the denominator polynomial be the same for every column of the Pade 
approximation. This forces the approximations of the unsteady aerodynamic 
forces that are dependent on the motion of a particular mode to be a function 
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of the variables that approximate the lag in development of the circulation for 
that mode and independent of the motion of the other modes. This results in a 
savings in the number of equations used to approximate the system because one 
set of equations are used to represent the lag in the development of the circu- 
lation per mode rather than per element as in the general Pade approximation. 


The problem is formulated as follows: 
mation in the frequency domain of the 
be a matrix of elements such that the 
polynomial of the Pade approximation; 
such that the jth diagonal element is 
approximation which best fits the jth 
Crete set of frequencies. Then, 

[Q ( jw) ] = [P(jw) ] [R(ju)) ] _1 


where 


let the matrix [Q(ju)) ] be the approxi- 
unsteady aerodynamic forces; let [PCjw) ] 
ith, jth element represents the numerator 
and let [R(jo>) ] be a diagonal matrix 
the denominator polynomial of the Pade 
column of the matrix fQ(ja)) ] at a dis- 


0 5) 


fR(jaj)] = (jw) N 



[I] 


* f w 

o=n 


0 6 ) 


N+2 

[P(jto) ] = ^ 
m=0 



rp m ] 


0 7) 


The details of the method of solution for finding the matrices [R^J and 
[P m ] of equations (16) and (17) are given in the appendix. The method is out- 
lined here. The problem is to find the matrices [P m ] and [R^] such that 
TQ(jw)] is a good approximation to TQ(jw)] subject to the constraints that 
the roots of the polynomials defined by [R] in equation (1 5) have negative 
real parts. This could be accomplished by minimizing 


J 



I 

n=1 


Nk 

^ i iQijdW 
n=i 


Qij (i k ii) ! 1 2 


0 8 ) 


while satisfying the stability constraints. In equation (18), J is the cost 
of the approximation, is the number of reduced frequencies at which 
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aerodynamic force data cure available, Qij (jkjl) is the ith,jth element of 
[Q(jk)] evaluated at the frequency k£, Qij(jk&) is the Pade approximation 
for the ith,jth element of [ Q ( j k) ] evaluated at the frequency k& , and 
|| || is a metric norm as defined in the appendix. By noting that the influ- 
ence of fp] and [r] on the cost is column dependent, by virture of [r] 
being diagonal, the problem can be reduced into n smaller problems by mini- 
mizing 


n N k 

jj = \ Yj E ^QijOkjl) - Qij(jk£)||2 ( 19 ) 

i=l JUl 

the function actually minimized in this paper is 
n Nk 

J j - ~ 2 ^ ^ liQij(jkJl) Rj(jkS,) “ Pij ( jk&) i I 2 (20) 

i=l &=1 


where Rj(jkjj,) is the jth element of the matrix [r] and Pij(jk£) is the 

ith,jth element of [p] evaluated at the reduced frequency k&. Minimizing 

equation (20) rather than equation (19) results in larger errors at the lower 
frequencies. It is shown subsequently to give good results. A numerical gra- 
dient procedure is used to find [R(ju)] so that equation (20) is minimized 

while constraining the roots of the polynomial to have negative real parts. 

At each gradient calculation, a linear least- squares estimator is used to cal- 
culate the elements of [P(jto)]. 

Equation (15) can be rewritten as 


[Q(jo>) ] 


[Aq] + (jo» 



[A]] + (jw) 2 



+ 




tR(jw) ] -1 


( 21 ) 


where [R(joo)] is equal to [R(ju) ] with the leading coefficient set to the 
identity matrix, or 
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[R(jw) ] 



CR ( j to) ] 


(22) 


Taking equation (10) and substituting CQ { jto) ] of equation (21) for [Q(jo>)J 
results in 

j([M] + £ C 2 [A 2 ]) (jw) 2 + ^[D] + ~ Vc [A-|])(jw) + ([K] + pV 2 [A 0 ]j 

N— "I 

1 V™' / c 

2° v2 Z W (; «*•>' 


m=0 


fR(jw) ] _1 / q(jw) = 0. 


The state-space formulation of this approximation to the aeroelastic model 
is derived by denoting the state variables as follows: 


"N 

x(t) t = F -1 fq(jw) ] 

x(t) 2 = F -1 [(jw) q(jo«] = x(t) -| 

x(t) 3 = P -1 [R -1 (jto) q(jo)) ] \ 

x(t) 4 - P-Htjw) R _1 (jw) q(jto) ] 

x(t) N+2 = F* 1 f{jto) N ” 1 R _1 (jto) q(jw) ] = x(t) N+1> 


(24) 


Using equation (24) , equation (23) can be rearranged as 


x(t)T = x ( t) 2 (25a) 

“([M] + ^ c 2 fA 2 ]) x ( t) 2 = ([K] + 2 V 2 [Aq]) x(t)T 

+ ( [D] + ^ Vc [A] ]) x (t) 2 
N-1 

1 / c \ro~N 

+ - pv 2 2 W (~) *(t)m*3 (25b) 

m=0 
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x(t ) 3 = x(t ) 4 


(25c) 


x(t) N+2 = x(t)-| 


N-l 

-I 

£=0 


/ c N 

{ — ) M x < fc H+3 


(25d) 


Equations (25) represent a set of first-order constant-coefficient differential 
equations which can be put into the form of 


X ( t) = A J X(t) = [c]" 1 [B] X(t) 


(26) 


The matrix [c] for N = 4 is 


[C] = 


[I] 


[m] + - c 2 [a 2 ] 
8 


0 

0 

0 

0 


[I] 

0 

0 

0 


0 

[l] 

£ 

0 


£ 

0 

[I] 

0 


£ 

0 

0 

[I] 


(27) 


and the matrix [b] for N = 4 is 


fB] 


[K] + - pV 2 [Aq] 
1 2 


[I] 


ru 


f D] + - VC [At 

4 


] • pv 2 [D 0 ](— 

2 'C 


2V\ < * 1 


2 “ Vc / 

0 
0 
0 

/2V\ 4 

- [R ° ! (r) 


( 2V\ 3 

PV 2 [D,] ( — ) 


[I] 

0 

0 


- t»ll (-) 


1 ( 2V\ 

- pv 2 [d 2 ] ( — ) 

2 'C f 


0 

[I] 

0 

/ 2V\ 2 

f*2l - 


1 ( 2V\ 

- PV 2 [D 3 ] ( — ) 

2 \ c ' 


- [r 3 ] 


(r) 


(28) 


/ PV^\ 

The dimensions of the matrix AlM,— I are a function of the number of modes 

used to represent the system n and the order of the Pade approximation used 
N. For clarity in the above development, the order of the approximation was 


14 



made the same for every mode shape in the system. This allowed the equations 
of motion to be derived with vector notation instead of using the individual 
scalar equation for each mode. In practice, the order of the approximation for 
each mode is adjusted until the error between the approximation and the aero- 
dynamic data is below a preset value. The number of first-order differential 
equations is equal to 


2n + 


2 

-i=l 


N-i 


(29) 


where N^ is the order of the Pade approximation used to approximate the 
unsteady aerodynamic forces for the ith mode. 

( PV2\ 

The eigenvalues of Al Mf ~ ~ ) include roots resulting from the unsteady- 

aerodynamic- forces approximation as well as the roots of the classical flutter 
equation. The value of the dynamic pressure pv^/2 that results in an oscil- 
latory eigenvalue with the real part of the eigenvalue equal to zero is the 
flutter point for the vehicle at the particular Mach number for which the 
unsteady aerodynamic data were calculated. Determining the flutter point for 
a range of Mach numbers determines the flutter boundary of the vehicle. 


RESULTS AND DISCUSSION 

The method developed in this paper is compared with the "Least-Squares 
Method" of references 3, 4, and 5 by application to the example of reference 5 
which is an aeroelastic semispan wind-tunnel model with an assumed plane of 
symmetry at the root. The wing geometry is given in figure 1 , and the general- 
ized masses and frequencies are presented in table I for the first 1 0 elastic 

modes. For the model used in this example, modes 1 , 2, 4, 5, and 6 were used 
to represent the wing. Structural damping was assumed to be negligible. The 
oscillatory aerodynamic forces were calculated using a doublet-lattice technique 
similar to that described in reference 8. In order to calculate the pressure 
distribution on an oscillating wing undergoing simple harmonic motion, the lift- 
ing surface is subdivided into an array of trapezoidal boxes arranged in strips 

parallel to the airstream as shown in figure 2. The lifting surface is then 

represented by a lattice of doublets located at the quarter chord of each box. 
The downwash condition is satisfied at the three-quarter chord of each box by 
equating it to the downwash resulting from the slope and deflection rate of each 
structural mode. The lifting surface was divided into 210 boxes arranged in 
30 strips spanwise with seven boxes chordwise. Oscillatory aerodynamic forces 
were calculated at six reduced frequencies (k = 0, 0.1, 0.3, 0.5, 0.7, and 0.9). 

The unsteady aerodynamic forces were approximated through the use of the 
method described in this present paper to calculate the matrix coefficients of 
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equation (21 ) and equation (22) . First- and second-order Pade approximations 
were used for all^five modes. The approximations (calculated in square meters) 
for $21 » $12' 311,3 $22 unsteady aerodynamic forces for N = 2 are 


Qll = (3.39 x 1 0 -8 ) ( joi) 2 + (1.18 x 1 0 -4 ) (jco) + (1.62 x 10" 2 ) ^ 

(2.03 x 1 0 3 ) (jco) (1 .59 x 1 0 2 ) ( jco) 

jco + 5.52 x 10 1 jco + 7.15 x 10 2 

$21 = (1.13 X 1 0“ 8 ) (jco) 2 - (2.28 X 10 _5 )(jco) + (2.56 x 1 0 -2 ) 

(5.25 x 10'3) (jco) (6.85 x 1 0“ 2 ) (jco) 

” + " " 

jco + 5.52 x 1 0 1 jco + 7.15 x 1 0 2 

$12 = (1.15 x 10-7) (jco) 2 - (1.51 x 10 _4 )(jco) - (1.69 x 1 0 _1 ) 

(3.42 x 10~ 2 ) (jco) (1.40 x 1 0 -2 ) ( jco) 
jco + 5.80 x 10^ jco + 4.26 x l 0 2 

$22 = - (8.47 x 1 0“8) (jco) 2 + (3.50 x 10 _4 )(jco) - (2.47 x 1 0 _1 ) 

(6.18 x 1 0~ 2 ) (jco) (8.64 x 1 0“ 2 ) (jco) 

+ + 

jco + 5.80 x 1 0 3 jco + 4.26 x l o 2 -s 


(30) 


The least-squares approximation for the same unsteady aerodynamic forces with 
the same a priori parameters as in reference 5 are 


Qn 


(3.81 x 1 0“8) (jco) 2 + (1.12 x 1 0“ 4 ) (jco) (1.63 x 1 0 -2 ) 


(5.37 x TO -3 ) (jco) 
jco + 8.50 x 1 0 1 


(8.49 x 10 -3 ) (jco) 
+ 

jco + 1.70 x i o 2 


(1 .23 x 1 0 -3 ) (jco) 
jco + 2.55 x 1 0 2 


(1.08 x 10 -2 ) (jco) 
jco + 3.40 x 1 0 2 


(31a) 
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Q 2 i = - (1.12 x 10" 9 ) ( jo >) 2 + (1.48 x l(T 6 )(jw) + (2.55 x 1 0~ 2 ) 

(1.56 x 1 0 -2 ) ( jw) (5.80 x 10“ 2 )(jw) 

— ' + * 

jw + 8.50 x TO 1 jw + 1.70 x 10 2 


( 1.47 x lO-’Hjw) (1 .37 x 10" 1 ) ( jco) 

+ (31b) 

joj + 2.55 x 10 2 jw + 3.40 x 10 2 


Ql2 58 “ (7.01 x 1 0 -9 ) ( ja>) 2 - (4.58 x 10 - 2 )(jw) - ( 1.68 x 1 0” 1 ) 
(7.45 x 1 0 -2 ) ( jco) (2.61 x 10 - 1 )(jw) 


jw + 8.50 x 10 1 jw + 1.70 x 10 2 
(6.41 x 10" 1 )(jw) (5.17 x 1 0 -1 ) ( jw) 

: — (31c) 

jw + 2.55 x 1 0 2 jw + 3.40 x 1 0 2 


Q22 = ~ (1 *50 * TO -7 ) (jw) 2 + ( 4.32 x 10" 4 )(jw) - (2.46 x 1 0 _1 ) 
(1.53 x 1 0 -1 ) ( jw) (3.93 x 10“ 1 )(jw) 


jw + 8.50 x 1 0 1 jw + 1.70 x 1 0 2 

( 6.88 x 1 0 - 1 ) ( jw) (3.69 x 1 0* 1 ) ( jw) 

+ — (31 d) 

jw + 2.55 x 1 0 2 jw + 3.40 x 1 0 2 


The fit, to the doublet-lattice data, achieved by the approximations described 
by equation (30) and equations (31) is presented in figure 3. For this case, 
figure 3 indicates that both the "Least-Squares Method" with four terms and the 
Pade approximation method presented in the present paper with two terms are 
accurate approximations for the Q 11f Q 21 , and Q 2 2 aerodynamic forces. The 
Pade approximation for the Qi 2 aerodynamic force does not do as well as the 
"Least-Squares Method." But, an important point to be noted is that the Pade 
method requires 33 percent fewer state equations than employed with the "Least- 
Squares Method." In addition, the a priori parameters required by the "Least- 
Squares Method" are not needed by the Pade approximation method. The frequen- 
cies and damping ratios predicted by the two methods are described in figure 4. 
Figure 4 also presents results of a first-order Pade approximation when used 
in the stability analysis. The figure indicates that the predicted damping 


17 



ratio for the first mode is not as accurate for the first-order Pade approxi- 
mation as for the second-order Pade approximation. Since the change in the 
frequency and the damping ratio with respect to a change in the dynamic pressure 
is predicted correctly near the flutter dynamic pressure, it may be possible to 
use a first-order Pade approximation in the design stage while adjusting the 
design parameters accordly to account for the error. If this could be done, 
a 50-percent reduction in the number of state equations needed to represent the 
system could be realized. As shown in figure 4, the second-order Pade approxi- 
mation estimates the dynamic pressure at which the system becomes unstable to 
about the same accuracy as does the "Least-Squares Method." 


CONCLUDING REMARKS 

An analytical method for approximating time-domain aerodynamic forces has 
been presented. The method is based on approximating the oscillatory aero- 
dynamic forces with Pade approximations in the reduced frequency domain. An 
approximate time-domain representation is then developed, assuming stable air- 
craft motion, through the use of the inverse Fourier transform. The analytical 
method is applied to an aeroelastic wind-tunnel model and showed good agreement 
with previously used analytical methods for predicting the flutter point and 
stability trends. Some of the important results of this work are as follows: 

1 . All the parameters of the aerodynamic model are calculated without a 
priori knowledge, unlike the currently formulated "Least-Squares Method." 

2. The resulting approximations are stable in the sense that a bounded 
input will produce a bounded generalized aerodynamic force output. 

3. The number of differential equations required to represent the aero- 
dynamic force to a given accuracy tends to be smaller than employed in certain 
existing techniques where the denominator coefficients are chosen a priori. 

The observed reduction in the set of system equations was 33 percent from 
that employed in the "Least-Squares Method." 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
September 5, 1980 
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APPENDIX 


APPROXIMATION OF UNSTEADY AERODYNAMIC FORCES 

The generalized aerodynamic forces can be thought of as a matrix of elements 
such that the ith,jth elements represent the generalized force that is applied 
to the ith structural mode due to a pressure distribution caused by the jth mode 
oscillating at a reduced frequency. In practice, the generalized aerodynamic 
force matrix is calculated for a number of specific Mach numbers and reduced 
frequencies. This appendix describes a technique which calculates an approxi- 
mation to the unsteady aerodynamic forces as a rational polynomial in reduced 
frequency. Details are only provided for the first-order Pade approximation 
(N = 1 ) . 

Consider equations (15), (16), and (17), with k being substituted for 

0)c/2V and N - 1 . Then 


CQ ( jk) ] = [p ( jk) ] [R(jJc) ]-"> 


(Al) 


where 


fo(jk) 3 = [r 0 3 + (jk) Cl ] 


(A2) 


and 


Tp(jk)] = [Pq] + (jk) fp n ] + (jk) 2 [p 2 ] + (jk) 3 [p 3 ] 


(A3) 


The matrices [Rq] , [Pq] , fPi ] , [P 2 ] t and Cp 3 ] must be found such that the 

following equation is made as small as possible and the roots of the polynomials 
of equation (A2) have negative real parts 


J 



I 

i=i 


I 

£=i 


I Q i j (jk&) - $ij(jkjl) II 2 


(A4) 


The function 


is the norm of the argument and is defined as 
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|x|| = \/(x,x) 


APPENDIX 


(A5) 


where the function ( , ) is the scalar product which has the following metric 
properties: 

(x,x) = 0 (equality holds if and only if x = 0 ) 

(y,x) = Conjugate (x,y) 

(x,cy) = c(x,y) (c = Constant) 

(x,yi+y 2 ) = (x,yi ) + (x,y 2 ) 


The minimization of equation (A4) can be separated into n smaller prob- 
lems by noting that the influence of [p] and [r] on the cost is column 
dependent since [r] is diagonal. Therefore, the problem is to find the jth 
column vectors [P-| r j], [P 2 f j], and [P 3 f j] of the matrices fPg] , 

[Pi ] , [p 2 ], and [P 3 ] and the jth element rg -i of the diagonal matrix 

[R 0 ] 

n Nk 

Jj = ^ ^ ||Qij(jkd) [r 0f j + (jk£)] - [P 0 ,j] - ( jk&) Oi,j] 

i=l JUl 

- (jkfc ) 2 fP 2 , j 1 - (jkg,) 3 I> 3 ,j ]|! 2 ( A6 ) 


The effect of minimizing equation (A 6 ) instead of equation (A4) is to cause the 
maximum relative error of the Pade approximation to be at the lower frequencies. 
For a given r 0,j' equation (A 6 ) can be rewritten as 




(A7) 


where 


APPENDIX 


R^[Qij(rQ^j — j^i ) 1 

Re [Qij( r O,j " 

Im[Qij (ro,j - jki ) ] 

Im [pij ( r 0, j - jk Nk 3 


A = 


lm[jki ] 


Im[jjk Nk "j 


Re Qjk] ) 0 

« 

Re H jk Nk)5 ° 

0 Im jj jkT ) 

* 

o im |jkn k )^ 


(A8) 


(A9) 


P 0,ij 

p l,ij 

p 2/ij 

P 3,ij 


(AT 0) 


where Pm, ij is th® ith, jth element of [P m ] . If the constraint that is 


| | Xj | | 2 = Minimum 


(All) 


imposed so that the requirement of the rank of the matrix A in equation (A9) 
does not have to equal the dimension of X j f the unique solution to equation (A7) 
exists as (ref. 9) 
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Xj = A + Bj (A12) 

where A^ is the pseudoinverse of A in equation (A9) . 

The minimization of equation (A6) is performed by a numerically constrained 
gradient procedure described in reference 10 for the variables of equation (A2) . 
The constraint that must be enforced is that the resulting approximation func- 
tion be stable. It is necessary and sufficient for the function to be stable 
if the roots of the functions described by equation (A2) have negative real 
parts. For the example presented in the appendix this requires that 


r 0,j 


< 0 


(A13) 


At each step in the parameter search for the minimum fo equation (A7) , equa- 
tion (A12) is used to solve for the variables in equation (A1 0) . The procedure 
for higher order approximations is similar. 
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TABLE I.- FREQUENCY 


Mode 


Natural frequency, 
Hz 


1 

2 

3 

4 

5 


5.233 
19.129 
20.906 
25.769 
46.110 
61 .234 
79.682 
86.030 
98.087 
118.150 


L 


GENERALIZED MASS 


Generalized mass 


kg 


3.678 
7.769 « 

7.044 
2.970 
4.714 
4.758 
5.156 
11.297 
7.558 
5.501 
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Imaginary part 




Least-squares approximation of 

reference 5 

o Present method, N = 1 
0 Present method, N = 2 


Frequency, Hz 


Damping ratio 



Figure 4.- Change in frequency and damping ratio with respect to change 

in dynamic pressure. 
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